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Abstract 


We introduce a RBFs mesheless method of lines that decomposes the 
interior and boundary centers to obtain the numerical solution of the time 
dependent PDEs. Then, the method is applied with an adaptive algorithm 
to obtain the numerical solution of one dimensional problems. We show that 
in the problems in which the solutions contain region with rapid variation, 
the adaptive RBFs methods are successful so that the PDE solution can be 
approximated well with a small number of basis functions. The method is 
described in detail, and computational experiments are performed for one- 
dimensional Burgers’ equations. 


Keywords: Method of Lines; Radial basis functions; Adaptive Method; 
Burgers’ Equations. 


1 Introduction 


The radial basis functions (RBFs) methods are one of the most attractive 
meshless methods. These methods are easy to implement, very suitable for 
problems in irregular geometries and the formulation for different dimensional 
problems are similar. Also, this method can be spectrally accurate [11]. A set 
of points called centers are needed to define the RBFs. Therefore, a RBF can 
be defined anywhere in a given domain, independently to the other RBFs. 
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Both the approximation quality and the stability of the RBFs interpolation 
depend on the positions of the centers set [9]. 

The condition number of RBFs collocation methods becomes large when 
the number of centers increases, while reducing the number of centers im- 
proves the conditioning [9, 13]. In order to obtain numerical solution with 
the minimal numbers of centers, we can use a set of adaptive nodes rather 
than uniform ones. Especially in problems whose solutions contain regions of 
rapid variation, adaptive methods are preferred over fixed grid methods, [17]. 
The goal of an adaptive method is to obtain a numerical solution such that 
the error is less than a prescribed accuracy but with the minimal number 
of grid points. By using adaptive methods, the computational grid should 
reflect the profile of the solution. Clearly, grids with finer spacing should be 
concentrated in regions, where high variations occur, and much coarser grids 
can be used in other regions. 

Some methods have been constructed to select centers of RBFs. In [6, 26], 
the power function is used to iteratively obtain an optimal set of nodes. In 
[25], an adaptive algorithm so-called residual sub-sampling is introduced such 
that nodes can be added or removed based on residuals evaluated at a finer 
set of nodes. Our goal is to move a fixed numbers of nodes in such a way that 
nodes move with time and concentrate in region of domain that the solution 
has rapid variations. To this goal in this paper, we use a simple adaptive 
nodes generation method that is used for finite difference computations [24] 
and RBFs method [23]. Also we introduce a RBF's meshless method of lines 
to solve time dependent PDE with adaptive centers. In this method, we 
divide centers to interior and boundary data centers and obtain the expansion 
coefficients of boundary centers as a function of interior ones. This gives an 
ODEs system that is only related to the expansion coefficients of the interior 
data centers instead of all data centers. Actually after approximation spatial 
derivatives of equation and boundary condition with RBFs, we have a system 
of differential algebraic equations (DAEs) [5]. By decomposing centers and 
replacing boundary coefficients as a function of interior ones we obtain a 
smaller system of ODEs. The resultant system of ODEs can be solved with 
a proper ODE solver. We use the function ode15s in Matlab for solving the 
resulting system of ordinary differential equations. 

In this paper, in order to combine the adaptive method and the RBFs 
method of lines, we start with a set of uniform centers, then the adaptive 
method is used to obtain new centers for initial condition. After obtaining 
the adaptive centers, the PDE is advanced for a small time step. The ode15s 
in Matlab is used for solving the resultant ODEs system. Then, the numer- 
ical solution of the PDE is used to obtain adaptive centers for next time. 
The procedure is repeated until the final time. We perform computational 
experiment for unsteady Burgers’ equations and demonstrate the benefits of 
adaptation in the numerical experiments. 

The rest of the paper is organized as follows. In Section 2 at first the 
RBFs method of lines is introduced, then adaptive method is extended for 
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time dependent PDEs. Numerical experiment are given in Section 3. Finally, 
the conclusion is given in Section 4. 


2 Meshless method of lines with adaptive RBFs 


In this section, we introduce a RBF's mesheless method of lines that decom- 
poses the interior and boundary centers to obtain the numerical solution of 
the time dependent PDEs. Then, apply the method with an adaptive algo- 
rithm to obtain the numerical solution of one dimensional Burgers’ equations. 


2.1 RBFs meshless method of lines 


There are two classes of RBF's, known as globally supported and locally 
supported [22, 16]. Globally supported RBFs are infinitely smoothed and 
contain a free parameter ¢€, called shape parameter. This parameter affects 
both accuracy of the solutions and conditioning of the collocation matrix. As 
e decreases, numerical solution of PDEs gets more accurate and the condition 
number of the resulting matrix gets larger. If the shape parameter becomes 
too small, the ill-conditioned matrix leads to numerical instabilities and loss 
of precision. Thus it is important to select a good values for e. There are 
some paper related to select an optimal value for RBFs shape parameter 
[21, 1, 14]. 

Generally a radial basis function is a function ¢;(a,¢) = ¢ (ella — x;|I2), 
which depends solely on the distance between x € R and a fixed center x; € 2. 
¢; : Rt > Ris a continuous function and ||-||2 represents the Euclidean norm. 
The multiquadrics (MQ) RBF proposed by Hardy [3], is one of the most 
used globally supported RBF's because of its spectral convergence property. 
In [4], Franke showed that the MQ RBF is one of the best methods among 
29 scattered data interpolation schemes. We here use MQ RBF defined as 
d(r,€) = f1+ (er). 

Let a set of N distinct centers {a;}%, is given in QU OQ, where 2 is a 
bounded domain in R. We assume that the arrangement of the centers is in 
such a way that the first N; centers and the last Ng centers lie in Q and OQ, 
respectively, N = N; + Ng. Consider the following time dependent PDE of 
the general form 


Ou(z, t) 
Ot 


—Lu(a#,t) = f(a,t), «red, Bu(ax,t) = g(x,t), x € OQ, 
(1) 


u(a#,0) = uo(x). (2) 


with the initial condition 
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£ and B are differential and boundary operators respectively. We approxi- 
mate the solution of equation (1) by 


Ny N 
uN (x) =S° ei(t)o(ja—aill) + S > ei(t)o (lla — 2il))- (3) 
i=1 i=Nr41 


Using collocation method to ensure that the approximation w’ (z) satisfies in 
equations (1), one obtains the following system of equations for the expansion 
coefficients 


AiiC + Ai2C2 = F + £L4(C1, C2), (4) 
0C, +0 C2 = Gt) — (A210 + Ao2C2), (5) 
where 
Ai i(i, 7) = O(|lvs—25||), i=1,...,Nr,9 =1,..., Nz, 
Ai o(i,j) = O(|lzs— all), i=1,--.,Nr,9 = Nri,.--,N, 
Aoi(i, 7) = Bo(||ai — ||), i= Nr4i,...,N,j =1,..., Nr, 
Ago(t, J) = Be(|lzi — 25[1), t= Niqi,.--, NG = Nrii,---,N, 
Lilie) = GiCr nC Cs), 
Ny N 
Lig(C1, C2) = > og (t)L¢ (lei — ayll) + SS (LO (lla - 25ll), 
j=l J=Nir4+1 


FT = [f (a1, t), ope ., f(ay,,t)], 
and 
G(t)* = [o(@Ny4154); soa ,g9(tn,t)]. 


Equations (4) and (5) are distinct from ODEs because the coefficient matrix 
of the CT = [C,, C4] is singular and are referred to as differential-algebraic 
equations (DAEs). DAEs differ in many ways from ordinary differential 
equations and there are some problems to be expected in solving these sys- 
tems. More information about differential-algebraic equations can be found 
in [10, 12]. In order to reach a system of ODEs, we obtain C2 and C> from 
equation (5) as follows: 


Co = Az> (G(t) — AgiCi), (6) 
Cy = AZ} (GH) - AaaCr). (7) 
Note that unlike the interpolation problem the invertibility of Az.2 may failed 


for some special centers arrangements. However, numerical experiments show 
that the cases of singularity for Kansa method is rare [19]. We substitute C2 
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and C2 into equation (4) to obtain a N; x Ny; nonlinear system of ordinary 
differential equation for C7; as follows: 


(Ali - A1,2.A3 52,1) O.=F+ Lo(C1) — Ai,2473G(t), (8) 


L4(C1)" = [Lig(Ci),---, £7 ¢(Cr)], 


Nr 
Lig(Ci) = S> cj (t)L¢ (llai — all) + s d;(t)L¢ (|x; — x5l)), 
j=l J=Ni41 


and d,(t) is jth component of the vector C2 = Az> (G(t) — AgiC}). 

After solving the reduced system using a proper ODE solver, its solu- 
tion vector C, is applied to obtain Cj and C’,, using the relations Cy = 
Az 5 (G(t) — A21C1) and C = [C,C]. This method can be used for high 
dimensional problems. In case of one dimensional problem, we have only two 
boundary nodes x; and ry. 


2.2 Adaptive method 


In this section, the proposed mesheless method of lines that decomposes 
the interior and boundary centers to reach a smaller system of equations is 
combined with an adaptive algorithm that is used for finite difference and 
RBFs computations [24, 23]. In this method, at first the arclength of the 
numerical solution is computed. Then, the total length is divided into (N—1) 
equal part and the projection of each part onto x-axis determines the position 
of adaptive centers. The selected nodes on z-axis are such that the variation 
of the solution is equi-distributed on each section. 

Suppose that the approximate solution and the centers are given at the 
time step t,. The adaptive method is generalized for RBFs and introduced 
in the following algorithm: 


1) Sy =0,$; = Sj + 4/(hB)?? + (w 


(ua? = a(t, 2"), A? = 8 — 24). 


This step compute the arclength of solution u at time step t”. 


2)5= 4 Wop k = 2,27 = 27,2) = 2h. 
In this step the total length is divided into (NV — 1) equal part. 


3) For j =2,---,N-1,A=(j-1)6. 


- while A> S, putk=k+1, 


ae (A—S,-1)h® ‘ 
- FP = Uy 1+ “gs *) Next j. 


These steps project each part on solution onto x-axis. 


54 A.R. Soheili, M. Arab Ameri and M. Barfeie 


The set %;,j = 2,--- ,(N — 1) are adaptive interior nodes and %1,%y are 
the boundary nodes which are fixed. In using adaptive centers in region with 
rapid variations, nodes are close to each other and hence a larger value of 
shape parameter is needed. In order to obtain results with a smaller shape 
parameter, the final set of centers are selected as .9%; + .12;. 


In solving PDE problems, at first we apply the above adaptive algorithm 
for the initial condition to obtain the adaptive centers at t = 0. Then, 
adaptive centers are used for the RBFs method of lines to advance the PDE 
for a small time step. Next, the approximate solution at this time is used 
to obtain the adaptive centers again. Note that in each step we need to 
interpolate u at the adaptive centers to obtain initial condition for next time. 
The procedure is repeated until approximate solution is obtained at the final 
time. 


3 Numerical experiments 


In this section, the proposed method is applied to obtain numerical solution 
of Burgers’ equation as follows: 


2 
Ou | wot = 1 _ 2 € (0,1), (9) 


where Re is the Reynolds number. Equation (9) has shock wave behavior 
when the coefficient of kinematic viscosity vy = 1/Re is small. Also, it is 
a useful model for many interesting physical problems such as modeling of 
fluid dynamics, turbulence, boundary layer behavior, shock wave formation, 
traffic flow and is an interesting test problem for establishing the efficiency 
of different methods [8, 20]. 

Example 1. We consider Equation (9) with the following exact solution [15] 


a+ w+ (u— a)jexp(n) 


1 + exp(n) 410) 


u(x,t) = 


where 7 = a.Re.(a — ut — B),a, uw and G are arbitrary constant. 
In this example a, 4 and £ are .4,.6 and .125 respectively. The boundary 
conditions are 


u(0,t)=1, u(l,t)=.2 t>0. (11) 


Initial condition is taken from the exact solution. In order to measure the 
error, root mean square error (rms) is computed at M evaluation nodes 2; 
as: 
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M 


» (u (z;) — u(zi))° 


M 


Table 1 shows the rms error at t = .2,.4,.6,.8 and t = 1 for Re = 100 and 
Re = 500. The results are computed for N = 50 adaptive centers. When Re 
increases, the gradient of solution become sharper and consequently a larger 
values of shape parameter is needed. In this example, the values of shape 
parameters for Re = 100 and Re = 500 are 50 and 150, respectively. 


The numerical solution in Example 1 at ¢ = .1,t = .5 and t = 1 for 
Re = 100 and Re = 500 are shown in Figures 1.a and 2.a respectively. 
Figures 1.b and 2.b show the corresponding nodes trajectories. Figures show 
that the nodes move with time and are concentrated in region with rapid 
variations. When Re increases, the gradient become sharper and the nodes 
are more concentrated in region with rapid variations. 


The numerical and exact solutions of Example 1 at t = 1 are plotted in 
Figure 3. In order to obtain numerical solution with a set of uniform centers a 
larger number of nodes is needed [23]. Figure 4 shows the numerical solutions 
and absolute errors for N = 50 uniform and adaptive centers. In the case 
of using uniform centers, the numerical solution with some oscillations is 
obtained for Re = 500 and « = 50 at t = .1. As Figure 5 shows in order 
to obtain an acceptebale solution at this time, we need to use more uniform 
nodes or a set of adaptive centers. 


Table 1: rms error values corresponding to Example 1 


Re rms error (t=.2) rms error (t=.4) rms error (t=.6) rms error (t=.8) _ rms error (t=1) 
100 2.121018¢-003 3.149610e-003 4.023452¢-003 4.826538¢-003 5.599954e-003 
500 1.485532e-003 1.840686e-003 2.517397e-003 4.152420e-003 6.497818¢e-003 


Example 2. We consider Burgers’ equation (9) with the initial condition 
u(a,0) = sin(ra), 
and the boundary conditions 


u(0,t) =u(1,t)=0, t>0. 


The exact solution for this example is given by [15] 


lund) Inv >>>, iA; sin(irx)exp(—i? 77) 
UL, = = . 
Ao + 0 j~, Ai cos(imx)exp(—i?1?vt) 


with the Fourier coefficients 
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° 
g 


Figure 1: The numerical solution and corresponding nodes trajectories for N = 50 and 
Re = 100 


0 0.2 0.4 0.6 0.8 1 


(a) (b) 
Figure 2: The numerical solution and corresponding nodes trajectories for N = 50 and 
Re = 500 
1 
Ao = | exp {—(2mv)7!(1 — cos(rx)} da, (13) 
0 
1 
Ae 2 | exp {—(2nv)"1(1 —cos(x))} cos(ime)dx, i>1. (14) 
0 


In this example, NV = 50 nodes are used. The computation are performed 
for a final time t = 3. The numerical solution at t = .01,t = .1,t=1,t=2 
and t = 3 for Re = 100 and Re = 500 are shown in Figures 6.a and 7.a 
respectively. Initial condition in Example 2 does not have rapid variation, 
but the variation of the solution increases with time. The variation increases 
until a time to less than ¢ = .75. After this time, the variation of the solution 
decreases. Nodes trajectories also have such behavior. The nodes trajectories 
are shown in Figures 6.b and 7.b. Nodes are moved with time and concen- 
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Figure 4: The numerical solutions (a) and absolute errors (b) in the case of using uniform 
and adaptive centers at t = .1 for Re = 500 


trated in region with rapid variations. For t < to, the variation increases and 
nodes are concentrated in region with rapid variation. For t > to, the varia- 
tion of the solution decreases with the time and hence the nodes trajectories 
diverge. 


Figure 8 shows the numerical solution, exact solution and the absolute 
error at t = 3 when Re = 100, « = 50 and N = 50. We can see that, the 
error of the proposed method method is as small as 1074. 


The numerical solution for Re = 500 are obtained for « = 110. In this 
case, obtaining numerical results with N = 50 and Re = 500 uniform centers 
is not possible as well. 
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Figure 5: The numerical and exact solution and absolute error at t = .1 when Re = 500 
for N = 100 uniform centers (top) and N = 50 adaptive centers (down) 


4 Conclusion 


The adaptive MQ RBF method of lines has been proposed for obtaining the 
numerical solution of Burgers’ equations. In the method of lines, centers 
in the domain were portioned into the interior and the boundary centers. 
By portioning centers and obtaining the expansion coefficients for boundary 
centers as a function of interior ones, the DAEs system was converted to a 
smaller ODEs system. The resulting ODE system was solved with ode15s 
in Matlab. Also, we have used a simple adaptive nodes generation method 
to enable the method for obtaining numerical solution of the problem with 
high gradient. In the adaptive method, the nodes moved with time and 
concentrated in region with rapid variation. When the gradient of solution 
increases the nodes become more closer in region with rapid variation. In 
this case numerical solution can be obtained with less number of centers in 
comparison with using uniform centers. 


Numerical experiments have been performed for one-dimensional Burgers’ 
equations. Numerical results show that the proposed adaptive method are 
preferred over fixed grid methods. For example, the adaptive method is able 
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u(xt) 


Figure 6: (a) The numerical solutions at different times for Re = 100 and (b) correspond- 
ing nodes trajectories in Example 2 


u(x,t) 


Figure 7: (a) The numerical solutions at different times for Re = 500 and (b) correspond- 
ing nodes trajectories in Example 2 


to solve Burgers’ equation for Re = 500 and N = 50 whereas the numerical 
solution could not be obtained for N = 50 uniform centers. 
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